arXiv: 1501.05701v3 [q-bio.PE] 9 Jun 2015 


Bayesian co-estimation of selling rate and locus-specific 
mutation rates for a partially selling population 

Benjamin D. Redelings, Seiji Kumagai, Liuyang Wang, 

Andrey Tatarenkov, Ann K. Sakai, Stephen G. Weller, 

Theresa M. Gulley John G. Avise, and Marcy K. Uyenoyama 


Abstract 

We present a Bayesian method for characterizing the mating system of populations reproducing through a 
mixture of self-fertilization and random outcrossing. Our method uses patterns of genetic variation across 
the genome as a basis for inference about pure hermaphroditism, androdioecy, and gynodioecy. We extend 
the standard coalescence model to accommodate these mating systems, accounting explicitly for multilocus 
identity disequilibrium, inbreeding depression, and variation in fertility among mating types. We incorporate 
the Ewens Sampling Formula (ESF) under the infinite-alleles model of mutation to obtain a novel expression 
for the likelihood of mating system parameters. Our Markov chain Monte Carlo (MCMC) algorithm assigns 
locus-specific mutation rates, drawn from a common mutation rate distribution that is itself estimated from 
the data using a Dirichlet Process Prior model. Among the parameters jointly inferred are the population¬ 
wide rate of self-fertilization, locus-specific mutation rates, and the number of generations since the most 
recent outcrossing event for each sampled individual. 


1. Introduction 


Inbreeding has pervasive consequences throughout the genome, affecting genealogical relationships between 
genes held at each locus within individuals and among multiple loci. This generation of genome-wide, 
multilocus disequilibria of various orders transforms the context in which evolution proceeds. Here, we 


address a simple form of inbreeding: a mixture of self-fertilization (selling) and random outcrossing (Clegg 
l9^|Ritland|[200^ . 


Various methods exist for the estimation of selfing rates from genetic data. Wright’s (19211 fundamental 


approach bases the estimation of selfing rates on the coefficient of inbreeding (Fjs), which reflects the 
departure from Hardy-Weinberg proportions of genotypes for a given set of allele frequencies. The maximum 


likelihood method of Enjalbert and David (20001 detects inbreeding from departures of multiple loci from 


Hardy-Weinberg proportions, estimating allele frequencies for each locus and accounting for correlations in 
heterozygosity among loci (identity disequilibrium, Cockerham and Weir||1968 1. David et al. (20071 extend 


the approach of Enjalbert and David (20001, basing the estimation of selfing rates on the distribution 


of heterozygotes across multiple, unlinked loci, while accommodating errors in scoring heterozygotes as 
homozygotes. A primary objective of Instruct ( Gao et ^|2007 l is the estimation of admixture. It extends 
the widely-used program structure (Pritchard et al. 20001, which bases the estimation of admixture on 


disequilibria of various forms, by accounting for disequilibria due to selfing. Progeny array methods (see 
Ritland|[2002 1, which base the estimation of selfing rates on the genetic analysis of progeny for which one or 
more parents are known, are particularly well-suited to plant populations. Wang et al. (20121 extend this 


approach to a random sample of individuals by reconstructing sibship relationships within the sample. 

Methods that base the estimation of inbreeding rates on the observed departure from random union 
of gametes require information on expected Hardy-Weinberg proportions. Population-wide frequencies of 
alleles observed in a sample at locus I ({pn}) can be estimated jointly in a maximum-likelihood framework 
(e.g., Hill et al. 19951 or integrated out as nuisance parameters in a Bayesian framework {e.g., Ayres and 
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Balding|[l998|. Similarly, locus-specific heterozygosity 




Pu 


( 1 ) 


can be obtained from observed allele frequencies ( Enjalbert and David|2000 | or estimated directly and jointly 
with the selfing rate ( David et aZ.||2007 1. 

In contrast, our Bayesian method for the analysis of partial self-fertilization derives from a coalescence 
model that accounts for genetic variation and uses the Ewens Sampling Formula (ESF, |Ewens||1972] |. Our 
approach replaces the estimation of allele frequencies or heterozygosity (Q by the estimation of a locus- 
specific mutation rate (0*) under the infinite-alleles model of mutation. We use a Dirichlet Process Prior 
(DPP) to determine the number of classes of mutation rates, the mutation rate for each class, and the class 
membership of each locus. We assign the DPP parameters in a conservative manner so that it creates a new 
mutational class only if sufficient evidence exists to justify doing so. Further, while other methods assume 
that the frequency in the population of an allelic class not observed in the sample is zero, the ESF provides 
the probability, under the infinite-alleles model of mutation, that the next-sampled gene represents a novel 


allele (see (22al). 


To estimate the probability that a random individual is uniparental (s*), we exploit identity disequilibrium 
(Cockerham and Weir 19681, the correlation in heterozygosity across loci. This association, even among 


unlinked loci, reflects that all loci within an individual share a history of inbreeding back to the most recent 
random outcrossing event. Conditional on the number of generations since this event, the genealogical 
histories of unlinked loci are independent. Our method infers the number of consecutive generations of self- 
fertilization in the immediate ancestry of each sampled diploid individual and the probability of coalescence 
during this period between the lineages at each locus. 

In inferring the full likelihood from the observed frequency spectrum of diploid genotypes at multiple 
unlinked loci, we determine the distributions of the allele frequency spectra ancestral to the sample at the 
most recent point at which all sampled gene lineages at each locus reside in separate individuals. At this 
point, the ESF provides the exact likelihood, obviating the need for further genealogical reconstruction. This 
approach permits computationally efficient analysis of samples comprising large numbers of individuals and 
large numbers of loci observed across the genome. 

Here, we address the estimation of inbreeding rates in populations undergoing pure hermaphroditism, an- 
drodioecy (hermaphrodites and males), or gynodioecy (hermaphrodites and females). Our method provides 
a means for the simultaneous inference of various aspects of the mating system, including the population 
proportions of sexual forms and levels of inbreeding depression. We apply our method to simulated data 
sets to demonstrate its accuracy in parameter estimation and in assessing uncertainty. Our application 


to microsatellite data from the androdioecious killifish Kryptolebias marmoratus (Mackiewicz et a/. 2006 


Tatarenkov et al]|2012| and to the gynodioecious Hawaiian endemic Schiedea salicaria (Wallace et a?.|2011| 


illustrates the formation of inferences about a number of biologically significant aspects, including measures 
of effective population size. 


Evolutionary model 


We describe our use of the Ewens Sampling Formula (ESF, Ewens|[T972 1 to determine likelihoods based on 
a sample of diploid multilocus genotypes. 

From a reduced sample, formed by subsampling a single gene from each locus from each diploid individ¬ 
ual, one could use the ESF to determine a likelihood function with a single parameter: the mutation rate. 


appropriately scaled to account for the acceleration of the coalescence rate caused by inbreeding (Nordborg 
and Donnelly||1997 Fu||1997 1. Observation of diploid genotypes provides information about another param¬ 
eter: the probability that a random individual is uniparental (uniparental proportion). We describe the 
dependence of these two composite parameters on the basic parameters of models of pure hermaphroditism, 
androdioecy, and gynodioecy. 
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Rates of coalescence and mutation 

Here, we describe the structure of the coalescence process shared by our models of pure hermaphroditism, 
androdioecy, and gynodioecy. 


Relative rates of coalescence and mutation: We represent the probability that a random individual is 
uniparental by s* and the probability that a pair of genes that reside in distinct individuals descend from the 
same parent in the immediately preceding generation by 1/A^*. These quantities determine the coalescence 
rate and the scaled mutation rate of the ESF. 

A pair of lineages residing in distinct individuals derive from a single parent (P) in the preceding gen¬ 
eration at rate 1/A^*. They descend from the same gene (immediate coalescence) or from distinct genes in 
that individual with equal probability. In the latter case, P is either uniparental (probability s*), implying 
descent once again of the lineages from a single individual in the preceding generation, or biparental, imply¬ 
ing descent from distinct individuals. Residence of a pair of lineages in a single individual rapidly resolves 
either to coalescence, with probability 


fc = 


S* 

2-s*’ 


( 2 ) 


or to residence in distinct individuals, with the complement probability. This expression is identical to the 
classical coefficient of identity (Wright 1921 Haldane 19241. The total rate of coalescence of lineages sampled 
from distinct individuals corresponds to 


(1 + / c )/2 


1 


N* N*{2-s*)' 

Our model assumes that coalescence and mutation occur on comparable time scales: 

lim ANu = 9 

N^OO 

u—^0 

lim N*/N = S, 

N—¥oo 

N *—>’00 


(3) 


(4) 


for u the rate of mutation under the infinite alleles model and N an arbitrary quantity that goes to infinity 
at a rate comparable to N* and 1/u. Here, S represents a scaled measure of effective population size (termed 


inbreeding effective size” by Crow and Denniston 19881, relative to a population comprising N reproductives. 


In large populations, switching of lineages between uniparental and biparental carriers occurs on the 
order of generations, virtually instantaneously relative to the rate at which lineages residing in distinct 


individuals coalesce (Nordborg and Donnelly 1997 Fu 19971. Our model assumes independence between 


the processes of coalescence and mutation and that these processes occur on a much longer time scale than 
random outcrossing: 

1-s* »u,l/iV*. (5) 

For m lineages, each residing in a distinct individual, the probability that the most recent event corresponds 
to mutation is 

mu 9* 


lim 


in which 


N^oo mu + (™) /[iV*(2 — s*)] 9* + m — 1’ 


N* 

9* = lim 2N*u{2 - s*) = lim ANu — (1 - s72) 

A^—foo AT—foo N 

li —>-0 u —>-0 


= 9{1 - s*/2)S, 

for 9 and S defined in ([^. In inbred populations, the single parameter of the FSF corresponds to 9*. 


(6) 
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Uniparental proportion and the rate of parent-sharing: In a population comprising Nh hermaphrodites, 
the rate of parent-sharing corresponds to 1/Nh, and the uniparental proportion (s*) corresponds to 


SH = 


ST -I- 1 — S 


(7a) 


for s the fraction of uniparental offspring at conception and r the rate of survival of uniparental relative to 
biparental offspring. For the pure-hermaphroditism model, we assign the arbitrary constant N in @ as Nh, 
implying 

Sh = 1. (7b) 

In androdioecious populations, comprising Nh reproducing hermaphrodites and Nm reproducing males 
(female-steriles), the uniparental proportion (s*) is identical to the case of pure hermaphroditism Q 


SA = 


(8a) 


ST -(- 1 — S 

A random gene derives from a male in the preceding generation with probability 

(1-sa)/2, 

and from a hermaphrodite with the complement probability. A pair of genes sampled from distinct individuals 
derive from the same parent (l/N*) with probability 


1 


[(l + s^)/2]2 , [{1-sa)I2Y 


Nh 


N„ 


(8b) 


In the absence of inbreeding (sa = 0), this expression agrees with the classical harmonic mean expression 
for effective population size (|Wright||1969||. For the androdioecy model, we assign the arbitrary constant in 
0 as the number of reproductives {Nh + A^m)) implying a scaled rate of coalescence corresponding to 


1 


Nh + Nm _ [(1 + s^)/2]^ [(1 — Sa)/‘^Y 


Na 


for 


1 Pm 
Nm 


Pm 


Pm — 


(8c) 


(9) 


Nh -f Nm 

the proportion of males among reproductive individuals. Relative effective number Sa S (0,1] takes its max¬ 
imum for populations in which the effective number Nj^, implied by the rate of parent sharing, corresponds 
to the total number of reproductives {Na = Nh + Nm)- At Sa = 1, the probability that a random gene 
derives from a male parent equals the proportion of males among reproductives: 

(1 'Sa)/2 — Pm- 

In gynodioecious populations, in which Nh hermaphrodites and Nf females (male-steriles) reproduce, the 
uniparental proportion (s*) corresponds to 


SG = 


rNhtt 


rNhtt + Nh{l — a) + Nfa ’ 


(10a) 


in which a represents the seed fertility of females relative to hermaphrodites and a the proportion of seeds 
of hermaphrodites set by self-pollen. A random gene derives from a female in the preceding generation with 
probability 

(1 - sg)F/2, 


F = 


Nfa 


Nh{l-a) + Nfa 


for 


(10b) 
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the proportion of biparental offspring that have a female parent. A pair of genes sampled from distinct 
individuals derive from the same parent (1/7V*) with probability 


1 [l-(l-SG)F/2p [{1-sg)F/2Y _ 

Nq Nh Nf 

We assign the arbitrary constant iV in Q as (Nh + Nf), implying a scaled rate of coalescence of 

1 Nh + Nf [1 - (1 - sg)F/2]^ [(1 - sg)F/2]^ 

So No 1 - Pf Pf 


(10c) 


(lOd) 


for 


Nf 

~ Nh + Nf 


( 11 ) 


the proportion of females among reproductive individuals. As for the androdioecy model, Sq S (0,1] achieves 
its maximum only if the proportion of females among reproductives equals the probability that a random 
gene derives from a female parent: 

{l-SG)F/2 = pf. 


Likelihood 

We here address the probability of a sample of diploid multilocus genotypes. 

Genealogical histories: For a sample comprising up to two alleles at each of L autosomal loci in n diploid 
individuals, we represent the observed genotypes by 


X={Xi,X2,...,Xi}, (12) 

in which X; denotes the set of genotypes observed at locus I, 

Xi = {XaX,2,...,Xz„}, (13) 

with 

Xife = {Xiki,Xik2) 

the genotype at locus I of individual k, with alleles Xiki and Xik 2 - 

To facilitate accounting for the shared recent history of genes borne by an individual in sample, we 
introduce latent variables 

T = {ri,T2,...,T„}, (14) 

for Tfc denoting the number of consecutive generations of selling in the immediate ancestry of the indi¬ 
vidual, and 

I = {hk}, (15) 

for Ilk indicating whether the lineages borne by the individual at locus I coalesce within the most recent 
Tk generations. Independent of other individuals, the number of consecutive generations of inbreeding in the 
ancestry of the k*^ individual is geometrically distributed: 

Tk ^ Geometric (s*), (16) 

with Tk = 0 signifying that individual k is the product of random outcrossing. Irrespective of whether 0, 1, 
or 2 of the genes at locus I in individual k are observed, Iik indicates whether the two genes at that locus in 
individual k coalesce during the Tk consecutive generations of inbreeding in its immediate ancestry: 


Ilk — 


0 if the two genes do not coalesce 
1 if the two genes coalesce. 
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a a P a a (3 7 * 



Ti = 2 T2 = 0 Ts = 3 r4 = 1 Ts = 0 

=0 72 = 0 73 = 1 74 = 1 75 = 0 

Figure 1 Following the history of the sample (X/) backwards in time until all ancestors of sampled genes 
reside in different individuals (Y;). Ovals represent individuals and dots represent genes. Blue lines indicate 
the parents of individuals, while red lines represent the ancestry of genes. Filled dots represent sampled 
genes for which the allelic class is observed (Greek letters) and their ancestral lineages. Open dots represent 
genes in the sample with unobserved allelic class (*). Grey dots represent other genes carried by ancestors 
of the sampled individuals. The relationship between the observed sample X; and the ancestral sample Y; 
is determined by the intervening coalescence events I;. T indicates the number of consecutive generations of 
selling for each sampled individual. 


Because the pair of lineages at any locus coalesce with probability ^ in each generation of selfing, 

Pr(/ifc = 0) = ^ = l-Pr(/ifc = l). (17) 

Figure [l] depicts the recent genealogical history at a locus I in 5 individuals. Individuals 2 and 5 are 
products of random outcrossing (T 2 = T 5 = 0 ), while the others derive from some positive number of 
consecutive generations of selfing in their immediate ancestry (Ti = = 3,T4 = 1). Both individuals 1 

and 3 are homozygotes (aa), with the lineages of individual 3 but not 1 coalescing more recently than the 
most recent outcrossing event {In = 0 , 7(3 = 1). As individual 2 is heterozygous (a/3), its lineages necessarily 
remain distinct since the most recent outcrossing event (J /2 = 0). One gene in each of individuals 4 and 5 
are unobserved (*), with the unobserved lineage in individual 4 but not 5 coalescing more recently than the 
most recent outcrossing event {I 14 = 1 , 7(5 = 0 ). 

In addition to the observed sample of diploid individuals, we consider the state of the sampled lineages 
at the most recent generation in which an outcrossing event has occurred in the ancestry of all n individuals. 
This point in the history of the sample occurs T generations into the past, for 

T = 1 + max Tk ■ 

k 

In Figure for example, 7 = 4, reflecting the most recent outcrossing event in the ancestry of individual 3. 
The ESF provides the probability of the allele frequency spectrum at this point. 

We represent the ordered list of allelic states of the lineages at T generations into the past by 

Y = {Yi,Y2,...,Yi}, (18) 


for Y( a list of ancestral genes in the same order as their descendants in X(. Each gene in Y( is the ancestor 
of either 1 or 2 genes at locus I from a particular individual in X( (131, depending on whether the lineages 
held by that individual coalesce during the consecutive generations of inbreeding in its immediate ancestry. 
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We represent the number of genes in Y/ by mi {n < mi < 2n). In Figure for example, X/ contains 10 
genes in 5 individuals, but Y; contains only 8 genes, with Yu the ancestor of only the first allele of Xji and 
Yi 5 the ancestor of both alleles of Xjs. 

We assume © that the initial phase of consecutive generations of selling is sufficiently short to ensure 
a negligible probability of mutation in any lineage at any locus and a negligible probability of coalescence 
between lineages held by distinct individuals more recently than T. Accordingly, the coalescence history I 


(151 completely determines the correspondence between genetic lineages in X (121 and Y (18l. 


Computing the likelihood: In principle, the likelihood of the observed data can be computed from the 
augmented likelihood by summation: 


Pr(X|0*, s*) = 5] ^ Pr(X, I, T|0*, s*). 


(19) 


for 


0 * = {91,e, 


2 ) ■ 


)} 


( 20 ) 


the list of scaled, locus-specific mutation rates, s* the population-wide uniparental proportion for the repro¬ 
ductive system under consideration {e.g., ^ for the pure hermaphroditism model), and T (14l and I (151 


the lists of latent variables representing the time since the most recent outcrossing event and whether the 
two lineages borne by a sampled individual coalesce during this period. Here we follow a common abuse of 
notation in using Pr(X) to denote Pr(X = x) for random variable X and realized value x. Summation (191 


is computationally expensive: the number of consecutive generations of inbreeding in the immediate ances¬ 
try of an individual (T^) has no upper limit (compare David et al. 2007| | and the number of combinations 
of coalescence states (Iik) across the L loci and n individuals increases exponentially (2^”) with the total 
number of assignments. We perform Markov chain Monte Carlo (MCMC) to avoid both these sums. 

To calculate the augmented likelihood, we begin by applying Bayes rule: 


Pr(X,I,T|0*,s*) = Pr(X,I|T,0*,s*)Pr(T|0* 


Because the times since the most recent outcrossing event T depend only on the uniparental proportion s* 
through (161, and not on the rates of mutation 0 *, 


Pr(T|0*,s*) = nPr(rfe|s*). 

fc=i 


Even though our model assumes the absence of physical linkage among any of the loci, the genetic data 
X and coalescence events I are not independent across loci because they depend on the times since the most 
recent outcrossing event T. Given T, however, the genetic data and coalescence events are independent 
across loci 

L 

Pr(X, I|T, 0*, s*) = n Pr(X,, I,|T, 91, s*). 

1^1 

Further, 


Pr{Xi,li\T,9:,s*) 


Pr{Xi\Ii,T,9:,s*)-PT{Ii\T,9t,s*) 

n 

Pr(Xj|I,,0;,s*)-nPr(/,fc|rfe). 

k=l 


This expression reflects that the times to the most recent outcrossing event T affect the observed genotypes 
X; only through the coalescence states 1; and that the coalescence states I; depend only on the times to the 
most recent outcrossing event T, through ([T7|. 
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To compute Pr(X;|Ii, 0*, s*), we incorporate latent variable Y; describing the states of lineages at 
the most recent point at which all occur in distinct individuals (Figure]^: 

Pt{Xi\1i,01s*) = Y.Pr{Xi,Yi\li,eis*) 

Y, 

= J2^iiXi\Yi,h,eis*)Pr{Yi\h,eis*) 

Y, 

= J2^iiXi\Yi,h)-PT{Yi\h,e:), (21a) 

Y, 

reflecting that the coalescence states 1/ establish the correspondence between the spectrum of genotypes in 
X; and the spectrum of alleles in Y; and that the distribution of Y;, given by the ESF, depends on the 
uniparental proportion s* only through the scaled mutation rate (|^. 

Given the sampled genotypes Xi and coalescence states I/, at most one ordered list of alleles Y; produces 
positive Pr(X;|Yi,I/) in ( |21a[ ). Coalescence of the lineages at locus I in any heterozygous individual {e.g., 
^ik = {Pi oi) with //fe = 1 in Figure]^ implies 

Pr(X,|Y,,I,)=0 


for all Y/. Any non-zero Pr(Xi|Y/,I;) precludes coalescence in any heterozygous individual and Y; must 
specify the observed alleles of X; in the order of observation, with either 1 (7;^ = 1) or 2 {In. = 0) instances 
of the allele for any homozygous individual {e.g., Xn. = {a, a)). For all cases with non-zero Pr(X;|Yi, I/), 


Pr(X,|Yj,I0 = l. 


Accordingly, expression (21a) reduces to 


p,{Xi\ii,eis*)= Y. P^{Yi\ii,0:), 

Y,:Pr(X,|Y,,I,)#0 


(21b) 


a sum with either 0 or 1 terms. Because all genes in Y; reside in distinct individuals, we obtain Pr(Y/|Ii, 0;*) 
from the Ewens Sampling Formula for a sample, of size 


mi 




k=l 


ordered in the sequence in which the genes are observed. 

To determine Pr(Yi|I;, 0^) in (21b), we use a fundamental property of the ESF (Ewens|1972 Karlin and 


McGreg^|l972 ): the probability that the next-sampled (z*^) gene represents a novel allele corresponds to 

0 * 

TTi = -—, (22a) 

i-i + 0 *' ^ ’ 

for 0 * defined in (1^, and the probability that it represents an additional copy of already-observed allele j is 


z — 1 


(22b) 


for ij the number of replicates of allele j in the sample at size (z — 1) (^^ ij = i — 1). Appendix K| presents 
a first-principles derivation of (22a). Expressions ( [2^ imply that for Y i the list of alleles at locu^ in order 
of observance, 

\ nil A — 1 1; 

(23) 


{ 0 *)Ki Y\Y{'mu-l)'- 

Pr(Yi|Ii,0r) = ^ 


nrJi(*-i+0r) ’ 

in which Ki denotes the total number of distinct allelic classes, mij the number of replicates of the allele 
in the sample, and m; = mij the number of lineages remaining at time T (Figure j^. 
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Missing data: Our method allows the allelic class of one or both genes at each locus to be missing. In 
Figure for example, the genotype of individual 4 is X /4 = (/3, *), indicating that the allelic class of the 
first gene is observed to be /3, but that of the second gene is unknown. 

A missing allelic specification in the sample of genotypes X; leads to a missing specification for the 
corresponding gene in Y/ unless the genetic lineage coalesces, in the interval between X; and Yp with a 
lineage ancestral to a gene for which the allelic type was observed. Figure illustrates such a coalescence 
event in the case of individual 4. In contrast, the lineages ancestral to the genes carried by individual 5 fail to 
coalescence more recently than their separation into distinct individuals, giving rise to a missing specification 
in Yz. 

The probability of Y; can be computed by simply summing over all possible values for each missing 
specification. Equivalently, those elements may simply be dropped from Y; before computing the probability 
via the ESF, the procedure implemented in our method. 


Bayesian inference framework 


Prior on mutation rates 


Ewens (19721 showed for the panmictic case that the number of distinct allelic classes observed at a locus 


{e.g., Ki in (231) provides a sufficient statistic for the estimation of the scaled mutation rate. Because each 
locus I provides relatively little information about the scaled mutation rate 6 ^ ([^, we assume that mutation 
rates across loci cluster in a finite number of groups. However, we do not know a priori the group assignment 
of loci or even the number of distinct rate classes among the observed loci. We make use of the Dirichlet 
process prior to estimate simultaneously the number of groups, the value of 9* for each group, and the 
assignment of loci to groups. 

The Dirichlet process comprises a base distribution, which here represents the distribution of the scaled 
mutation rate 9* across groups, and a concentration parameter a, which controls the probability that each 
successive locus forms a new group. We assign 0.1 to a of the Dirichlet process, and place a gamma 
distribution (r(Q; = 0.25,/? = 2)) on the mean scaled mutation rate for each group. As this prior has a high 
variance relative to the mean (0.5), it is relatively uninformative about 9*. 


Model-specific parameters 

Derivations presented in the preceding section indicate that the probability of a sample of diploid genotypes 
under the infinite alleles model depends on only the uniparental proportion s* and the scaled mutation rates 
0* (20 1 across loci. These composite parameters are determined by the set of basic demographic parameters 


tF associated with each model of reproduction under consideration. As the genotypic data provide equal 
support to any combination of basic parameters that implies the same values of s* and 0*, the full set 
of basic parameters for any model are in general non-identifiable using the observed genotype frequency 
spectrum alone. 

Even so, our MCMC implementation updates the basic parameters directly, with likelihoods determined 
from the implied values of s* and 0*. This feature facilitates the incorporation of information in addition to 
the genotypic data that can contribute to the estimation of the basic parameters under a particular model 
or assessment of alternative models. We have 


Pr(X, 0*, = Pr(X|0*, ’S') • Pr(0*) • Pr(’®') 

= Pr(X|0*,s*(’®')) • Pr(0*) • Pr(’®'), 


(24) 


for X the genotypic data and s*(’®') the uniparental proportion determined by ’S' for the model under 
consideration. To determine the marginal distribution of 9i (|^ for each locus I, we use (|^, incorporating 
the distributions of s*(’®') and S' (S'), the scaling factor defined in (|^: 


0i = 


S{l-s*/2)' 
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For the pure hermaphroditism model 0, = {s, r}, where s is the proportion of conceptions through 

selling, and r is the relative viability of uniparental offspring. We propose uniform priors for s and r: 


s ^ Uniform(0,1) 
T ^ Uniform(0, 1). 


(25) 


For the androdioecy model 0, we propose uniform priors for each basic parameter in fl/ = {s,T,Pm}' 

s ~ Uniform(0,1) 

T ~ Uniform(0,1) (26) 

Pm ~ Uniform(0,1). 


For the gynodioecy model (10 1 , ' 5 ' = {a,T,pf,a}, including a the proportion of egg cells produced by 
hermaphrodites fertilized by selling, pf ( 11 ) the proportion of females (male-steriles) among reproductives, 
and a the fertility of females relative to hermaphrodites. We propose the uniform priors 


a ^ Uniform(0,1) 
T ^ Uniform(0,1) 
Pf ^ Uniform(0,1) 
l/cr ~ Uniform(0,1). 


(27) 


Assessment of accuracy and coverage using simulated data 


We developed a forward-in-time simulator (https;//github. com/skumagai/selfingsim) that tracks mul¬ 
tiple neutral loci with locus-specific scaled mutation rates (©) in a population comprising N reproducing 
hermaphrodites of which a proportion s* are of uniparental origin. We used this simulator to generate data 
under two sampling regimes: large (L = 32 loci in each of n = 70 diploid individuals) and small (L = 6 loci 


in each of n = 10 diploid individuals). We applied our Bayesian method and RMES (David et al. 2007) to 


simulated data sets. A description of the procedures used to assess the accuracy and coverage properties of 
the three methods is included in the Supplementary Online Material. 

In addition, we determine the uniparental proportion (s*) inferred from the departure from Hardy- 
Weinberg expectation (F/g, Wright|[l969 1 alone. Our F/g-based estimate entails setting the observed value 
of Fjs equal to its classical expectation s */(2 — s*) ( [Wright |[l 92 1[ |Haldane||1924[ ) and solving for s*: 


‘^Fjs 
1 + Fjs 


(28) 


In accommodating multiple loci, this estimate incorporates a multilocus estimate for Fjs (Appendix]^ but, 
unlike those generated by our Bayesian method and RMES, does not use identity disequilibrium across loci 
within individuals to infer the number of generations since the most recent outcross event in their ancestry. 


As our primary purpose in examining the F/s-based estimate (28) is to provide a baseline for the results of 


those likelihood-based methods, we have not attempted to develop an index of error or uncertainty for it. 


Accuracy 

To assess relative accuracy of estimates of the uniparental proportion s*, we determine the bias and root- 
mean-squared error of the three methods by averaging over 10 "^ data sets ( 10 ^ independent samples from each 
of 10^ independent simulations for each assigned s*). In contrast with the point estimates of s* produced 
by RMES, our Bayesian method generates a posterior distribution. To facilitate comparison, we reduce our 
estimate to a single value, the median of the posterior distribution of s*, with the caveat that the mode and 
mean may show different qualitative behavior (see Supplementary Online Material). 

Figure [^indicates that both RMES and our method show positive bias upon application to data sets for 
which the true uniparental proportion s* is close to zero and negative bias for s* close to unity. This trend 
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Figure 2 Errors for the full likelihood (posterior median), RMES, and F/g-based (28 1 methods for a large 
simulated sample (n = 70 individuals, L = 32 loci). In the legend, rms indicates the root-mean-squared 
error and bias the average deviation. Averages are taken across simulated data sets at each true value of s*. 


reflects that both methods yield estimates of s* constrained to lie between 0 and 1. In contrast, the Fjs- 
based estimate (281 underestimates s* throughout the range, even near s* = 0 (Fjs is not constrained to be 
positive). Our method has a bias near 0 that is substantially larger than the bias of RMES, and an error that 
is slightly larger. A major contributor to this trend is that our Bayesian estimate is represented by only the 
median of the posterior distribution of the uniparental proportion s*. Figure [^indicates that for data sets 
generated under a true value of s* of 0 (full random outcrossing), the posterior distribution for s* has greater 
mass near 0. Further, as the posterior mode does not display large bias near 0 (Figure [Slj), we conclude that 
the bias shown by the median (Figure merely represents uncertainty in the posterior distribution for s* 
and not any preference for incorrect values. We note that our method assumes that the data are derived 
from a population reproducing through a mixture of self-fertilization and random outcrossing. Assessment 
of a model of complete random mating (s* = 0) against the present model (s* > 0) might be conducted 
through the Bayes factor. 

Except in cases in which the true s* is very close to 0, the error for RMES exceeds the error for our method 
under both sampling regimes (Figure]^. RMES differs from the other two methods in the steep rise in both 
bias and rms error for high values of s*, with the change point occurring at lower values of the uniparental 
proportion s* for the small sampling regime (n = 10, L = 6). A likely contributing factor to the increased 
error shown by RMES under high values of s* is its default assumption that the number of generations in 
the ancestry of any individual does not exceed 20. Violations of this assumption arise more often under 
high values of s*, possibly promoting underestimation of the uniparental proportion. Further, RMES discards 
data at loci at which no heterozygotes are observed, and terminates analysis altogether if the number of loci 
drops below 2. RMES treats all loci with zero heterozygosity 0 as uninformative, even if multiple alleles 
are observed. In contrast, our full likelihood method uses data from all loci, with polymorphic loci in the 
absence of heterozygotes providing strong evidence of high rates of selling (rather than low rates of mutation). 
Under the large sampling regime (n = 70, L = 32), RMES discards on average 50% of the loci for true s* 
values exceeding 0.94, with less than 10% of data sets unanalyzable (fewer than 2 informative loci) even at 
s* = 0.99 (Figure]^. Under the n = 10, L = 6 regime, RMES discards on average 50% of loci for true s* 
values exceeding 0.85, with about 50% of data sets unanalyzable under s* > 0.94. 

The error for the F/g-based estimate (281 also exceeds the error for our method. It is largest near s* = 0 
and vanishes as s* approaches 1, a pattern distinct from RMES (Figure]^. 
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Average posterior density when s* = 0 



Figure 3 Average posterior density of the uniparental proportion (s*) inferred from simulated data generated 
under the large sample regime (n = 70, L = 32) with a true value of s* = 0. The average was taken across 
posterior densities for 100 data sets. 
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Figure 4 Fraction of loci and data sets that are ignored by RMES. 


Coverage 

We determine the fraction of data sets for which the confidence interval (Cl) generated by RMES and the 
Bayesian credible interval (BCI) generated by our method contains the true value of the uniparental pro¬ 
portion s*. This measure of coverage is a frequentist notion, as it treats each true value of s* separately. 
A 95% Cl should contain the truth 95% of the time for each specific value of s*. However, a 95% BCI is 
not expected to have 95% coverage at each value of s*, but rather 95% coverage averaged over values of s* 
sampled from the prior. Of the various ways to determine a BCI for a given posterior distribution, we choose 
to report the highest posterior density BCI (rather than the central BCI, for example). 

Figure indicates that coverage of the 95% CIs produced by RMES are consistently lower than 95% 
across all true s* values under the large sampling regime {n = 70 L = 32). Coverage appears to decline 
as s* increases, dropping from 86% for s* = 0.1 to 64% for s* = 0.99. In contrast, the 95% BCIs have 
slightly greater than 95% frequentist coverage for each value of s*, except for s* values very close to the 
extremes (0 and 1). Under very high rates of inbreeding (s* « 1), an assumption ([^ of our underlying model 
(random outcrossing occurs on a time scale much shorter than the time scales of mutation and coalescence) 
is likely violated. We observed similar behavior under nominal coverage levels ranging from 0.5 to 0.99 
(Supplementary Material). 
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Figure 5 Frequentist coverage at each level of s* for 95% intervals from RMES and the method based on the 
full likelihood under the large sampling regime (n = 70,L = 32). RMES intervals are 95% confidence intervals 
computed via profile likelihood. Full likelihood intervals are 95% highest posterior density Bayesian credible 
intervals. 


Number of consecutive generations of selfing 

In order to check the accuracy of our reconstructed generations of selfing, we examine the posterior distri¬ 
butions of selfing times {T^} for s* = 0.5 under the large sampling regime (n = 70, L = 32). We average 
posterior distributions for selfing times across 100 simulated data sets, and across individuals k = 1... 70 



Figure 6 Exact distribution of selfing times under s* = 0.5 compared to the posterior distribution averaged 
across individuals and across data sets. 

within each simulated data set. We then compare these averages based on the simulated data with the 
exact distribution of selfing times across individuals (Figure]^. The pooled posterior distribution closely 
matches the exact distribution. This simple check suggests that our method correctly infers the true posterior 
distribution of selfing times for each sampled individual. 
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Analysis of microsatellite data from natural populations 

Androdioecious vertebrate 


Our analysis of data from the androdioecious killifish Kryptolebias marmoratus (Mackiewicz et al. 2006 


Tatarenkov et al. 20121 incorporates genotypes from 32 microsatellite loci as well as information on the 


observed fraction of males. Our method simultaneously estimates the proportion of males in the population 
{Pm) together with rates of locus-specific mutation (9*) and the uniparental proportion (s^). We apply the 
method to two populations, which show highly divergent rates of inbreeding. 


Parameter estimation: Our androdioecy model (251 comprises 3 basic parameters, including the fraction 
of males among reproductives (pm) and the relative viability of uniparental offspring (r). Our analysis 
incorporates the observation of rim males among ritotai zygotes directly into the likelihood expression: 


FT{X,l,T,nm\s* ,Q*,Pm,ntotal) = Pr(X,I,T|s*, 0*) • PT{nm\Pm,ntotal), 


in which 

~ Binomial(ntotoz,pm), (29) 

reflecting that s* and 0* are sufficient to account for X, I, and T, and also independent of rim, ntotai, and 
Pm- 

In the absence of direct information regarding the existence or intensity of inbreeding depression, we 
impose the constraint r = 1 to permit estimation of the uniparental proportion sa under a uniform prior: 

s* ^ Uniform(0,1). 


Low outcrossing rate: We applied our method to the BP data set described by Tatarenkov et al. (20121. 
This data set comprises a total of 70 individuals, collected in 2007, 2010, and 2011 from the Big Pine location 
on the Florida Keys. 


Tatarenkov et al. (2012| report 21 males among the 201 individuals collected from various locations 
in the Florida Keys during this period, consistent with other estimates of about 1% {e.g., Turner et al. 


19921. Based on the long-term experience of the Tatarenkov-Avise laboratory with this species, we assumed 
observation of n,„ = 20 males out of ritotai = 2000 individuals in (29l. We estimate that the fraction of 
males in the population (pm) has a posterior median of 0.01 with a 95% Bayesian Credible Interval (BCI) 
of (0.0062,0.015). 

Our estimates of mutation rates {9*) indicate substantial variation among loci, with the median ranging 
over an order of magnitude (ca. 0.5-5.0) (Figure [S^ Supplementary Material). The distribution of mutation 
rates across loci appears to be multimodal, with many loci having a relatively low rate and some having 
larger rates. 

Figure [^shows the posterior distribution of uniparental proportion s^, with a median of 0.95 and a 95% 
BCI of (0.93,0.97). This estimate is somewhat lower than F/s-based estimate ( [28| of 0.97, and slightly 
higher than the RMES estimate of 0.94, which has a 95% Confidence Interval (Cl) of (0.91,0.96). We note 
that RMES discarded from the analysis 9 loci (out of 32) which showed no heterozygosity, even though 7 of 
the 9 were polymorphic in the sample. 

Our method estimates the latent variables {Ti, T 2 ,..., T„} ( |14[ ), representing the number of generations 
since the most recent outcross event in the ancestry of each individual (Figure S5). Figure shows the 
empirical distribution of the time since outcrossing across individuals, averaged over posterior uncertainty, 
indicating a complete absence of biparental individuals (0 generations of selling). Because we expect that a 
sample of size 70 would include at least some biparental individuals under the inferred uniparental proportion 
(sA « 0.95), this finding suggests that any biparental individuals in the sample show lower heterozygosity 
than expected from the observed level of genetic variation. This deficiency suggests that an extended model 
that accommodates biparental inbreeding or population subdivision may account for the data better than 
the present model, which allows only selling and random outcrossing. 
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Figure 7 Posterior distribution of the uniparental proportion sa for the BP population. The median is 
indicated by a black dot, with a red bar for the 95% BCI and an orange bar for the 50% BCI. 
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Figure 8 Empirical distribution of number of generations since the most recent outcross event (T) across 
individuals for the K. marmoratus (BP population), averaged across posterior samples. The right panel 
is constructed by zooming in on the panel on the left. “Expected” probabilities represent the proportion 
of individuals with the indicated number of selling generations expected under the estimated uniparental 
proportion sa- “Inferred” probabilities represent proportions inferred across individuals in the sample. The 
first inferred bar with positive probability corresponds to T = 1. 


Higher outcrossing rate: We apply the three methods to the sample collected in 2005 from Twin Cays, 
Belize (TC05: Mackiewicz et al. 20061. This data set departs sharply from that of the BP population, 
showing considerably higher incidence of males and levels of polymorphism and heterozygosity. 

We incorporate the observation of 19 males among the 112 individuals collected from Belize in 2005 
( Mackiewicz et al|2006 1 into the likelihood (see (291). Our estimate of the fraction of males in the population 
(Pm) has a posterior median of 0.17 with a 95% BCI of (0.11,0.25). 

Figure S 6 (Supplementary Material) indicates that the posterior medians of the locus-specific mutation 
rates range over a wide range (ca. 0.5-23). Two loci appear to exhibit a mutation rates substantially higher 
than other loci, both of which appear to have high rates in the BP population as well (Figure [S4|. 

All three methods confirm the inference of |Mackiewicz et^. (20061 of much lower inbreeding in the TC 
population relative to the BP population. Our posterior distribution of uniparental proportion has a 
median and 95% BCI of 0.35 (0.25, 0.45) (Figure]^. The median again lies between the F/ 5 -based estimate 
(281 of 0.39 and the RMES estimate of 0.33, with its 95% Cl of (0.30,0.36). In this case, RMES excluded from 
the analysis only a single locus, which was monomorphic in the sample. 

Figure shows the inferred distribution of the number of generations since the most recent outcross 
event (T) across individuals, averaged over posterior uncertainty. In contrast to the BP population, the 
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Figure 9 Posterior distribution of the uniparental proportion for the TC population. Also shown are the 
95% BCI (red), 50% BCI (orange), and median (black dot). 




Generations 

Figure 10 Empirical distribution of selfing times T across individuals, for K. marmoratus (Population TC). 
The histogram is averaged across posterior samples. 


distribution of selfing time in the TC population appears to conform to the distribution expected under the 
inferred uniparental proportion (sa)> including a high fraction of biparental individuals (T^ = 0). Figure 
S7 (Supplementary Material) presents the posterior distribution of the number of consecutive generations of 
selfing in the immediate ancestry of each individual. 


Gynodioecious plant 

We next examine data from Schiedea salicaria, a gynodioecious member of the carnation family endemic 
to the Hawaiiian islands. We analyzed genotypes at 9 microsatellite loci from 25 S. salicaria individuals 


collected from west Maui and identified by Wallace et at. (20111 as non-hybrids. 


Parameter estimation: Our gynodioecy model (271 comprises 4 basic parameters, including the relative 


seed set of females (a) and the relative viability of uniparental offspring (r). Our analysis of microsatellite 


data from the gynodioecious Hawaiian endemic Schiedea salicaria ( 

Wallace et al. 

20111 constrained the 

relative seed set of females to unity (cr = 1), consistent with empirical results ( 

Weller and Sakai 

2005 

1 . In 


addition, we use results of experimental studies of inbreeding depression to develop an informative prior 
distribution for r: 

r~Beta(2,8), (30) 


the mean of which (0.2) is consistent with the results of greenhouse experiments reported by Sakai et al. 


(19891. 
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Campbell et al. (20101 reported a 12% proportion of females {uf = 27 females among ritotai = 221 
individuals). As in the case of androdioecy (|2^, we model this information by 


Uf ~ Binomial(ntoto;,p/), 


(31) 


obtaining estimates from the extended likelihood function corresponding to the product of PT{nf\ntotai,Pf) 
and the likelihood of the genetic data. We retain a uniform prior for the proportion of seeds of hermaphrodite 
set by self-pollen (a). 


Results: Figure SIO (Supplementary Material) presents posterior distributions of the basic parameters 
of the gynodioecy model (lOl. Our estimate of the uniparental proportion sq (median 0.247, 95% BCI 

of Sc = 0.33. Although RMES excluded 


(.0791, 0.444)) is substantially lower than the F/g-based estimate 
none of the loci, it gives an estimate of sq = 0, with a 95% Cl of (0, 0.15). 

Unlike the K. marmoratus data sets, the S. salicaria data set does not appear to provide substantial 
evidence for large differences in locus-specific mutation rates across loci: Figure [S^ (Supplementary Material) 
shows similar posterior medians for across loci. 

Figure mi presents the inferred distribution of the number of generations since the most recent outcross 



Figure 11 Empirical distribution of selfing times T across individuals, for S. salicaria. The histogram is 
averaged across posterior samples. 

event T across individuals, averaged over posterior uncertainty. In contrast with the analysis of the K. 
marmoratus BP population (Figure]^, the distribution appears to be consistent with the inferred uniparental 
proportion sc- Figure (Supplementary Material) presents the posterior distribution of the number of 
consecutive generations of selfing in the immediate ancestry of each individual. 

Table [^presents posterior medians and 95% BCIs for the proportion of uniparentals among reproductives 
(s*), the proportion of seeds set by hermaphrodites by self-pollen (a), the viability of uniparental offspring 
relative to biparental offspring (r), the proportion of females among reproductives (p/), and the probability 
that a random gene derives from a female parent ((1 — sc)F/2). Comparison of the first (YYY) and 
fifth (NYY) rows indicates that inclusion of the genetic data more than doubles the posterior median of s* 
(from 0.112 to 0.247) and shrinks the credible interval. Comparison of the first (YYY) and third (YNY) rows 
indicates that counts of females and hermaphrodites greatly reduce the posterior median of p/ and accordingly 
change the proportional contribution of females to the gene pool ((1 — sq)F/2). The bottom row of the table 
(NNN), showing a prior estimate for composite parameter s* of 0.0844(0.000797,0.643), illustrates that its 
induced prior distribution departs from uniform on (0,1), even though both of its components (a and r) 
have uniform priors. 










Table 1 Parameter estimates for different amounts of data. Estimates are given by a posterior median and a 95% BCI. 
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Discussion 

We introduce a model-based Bayesian method for the inference of the rate of self-fertilization and other 
aspects of a mixed mating system. In anticipation of large (even genome-scale) numbers of loci, it uses 
the Ewens Sampling Formula (ESF) to determine likelihoods in a computationally efficient manner from 
frequency spectra of genotypes observed at multiple unlinked sites throughout the genome. Our MCMC 
sampler explicitly incorporates the full set of parameters for each iconic mating system considered here 
(pure hermaphroditism, androdioecy, and gynodioecy), permitting insight into various components of the 
evolutionary process, including effective population size relative to the number of reproductives. 


Assessment of the new approach 


Accuracy: Enjalbert and David| (|2000| and [David et al. (20071 base estimates of selling rate on the distribu¬ 


tion of numbers of heterozygous loci. Both methods strip genotype information from the data, distinguishing 
between only homozygotes and heterozygotes, irrespective of the alleles involved. Loci lacking heterozygotes 
altogether (even if polymorphic) are removed from the analysis as uninformative about the magnitude of 
departure from Hardy-Weinberg proportions (Figure]^. As the observation of polymorphic loci with low 
heterozygosity provides strong evidence of inbreeding, exclusion of such loci by RMES ( David et a/. [2007 1 may 
contribute to its loss of accuracy for high rates of selling (Figure]^. 

Our method derives information from all loci. Like most coalescence-based models, it accounts for the 
level of variation as well as the way in which variation is partitioned within the sample. Even a locus 
monomorphic within a sample provides information about the age of the most recent common ancestor of 
the observed sequences, a property that was not widely appreciated prior to analyses of the absence of 
variation in a sample of human Y chromosomes ( Dorit et a7||1995l|Fu and Li||1996 l. 

Estimates of the rate of inbreeding produced by our method appear to show greater accuracy than RMES 
and the F/s-based method (281 over much of the parameter range (Figure]^. The increased error exhibited 
under very high rates of inbreeding (s* « 1) may reflect violation of our assumption ([^ that random 
outcrossing occurs on a much shorter time scale than mutation and coalescence. Even though our method 
assumes that the rate of inbreeding lies in (0,1), the posterior distribution for data generated under random 
outcrossing (s* = 0) does indicate greater confidence in low rates of inbreeding (Figure]^. 

Both RMES and our method invoke independence of genealogical histories of unlinked loci, conditional on 
the time since the most recent outcrossing event. RMES seeks to approximate the likelihood by summing over 
the distribution of time since the most recent outcross event, but truncates the infinite sum at 20 generations. 
The increased error exhibited by RMES under high rates of inbreeding may reflect that the likelihood has a 
substantial mass beyond the truncation point in such cases. Our method explicitly estimates the latent 


variable of time since the most recent outcross for each individual (141. This quantity ranges over the non¬ 


negative integers, but values assigned to individuals are explored by the MCMC according to their effects 
on the likelihood. 


Frequentist coverage properties: Bayesian approaches afford a direct means of assessing confidence in pa¬ 
rameter estimates, and our simulation studies suggest that the Bayesian Credible Intervals (BCIs) generated 
by our method have relatively good frequentist coverage properties as well (Figure S3). The Confidence 
Intervals (CIs) reported by the maximum-likelihood method RMES (David et al. ]2007| appear to perform less 
well (Figure [5|. Although David et al. (2007) describe RMES as determining CIs via the profile likelihood 
method (see jKreutz et a7|2013 |, RMES holds constant parameters other than the uniparental proportion (s*) 
instead of reoptimizing them to maximize the likelihood as s* varies. The result is therefore not a true profile 
likelihood, which may explain the poor coverage properties of the CIs that RMES provides. 


Model fit: Bayesian approaches also afford insight into the suitability of the underlying model. Our method 
provides estimates of the number of generations since the most recent outcross event in the immediate 
ancestry of each individual (T). We can pool such estimates of selling times to obtain an empirical distribution 
of the number of selling generations, a procedure particularly useful for samples containing observation of the 
genotype of many individuals. Under the assumption of a single population-wide rate of self-fertilization, we 
expect selling time to have a geometric distribution with parameter corresponding to the estimated selling 
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rate. Empirical distributions of the estimated number of generations since the last outcross appear consistent 
with this expectation for the data sets derived from the TC population of K. marmoratus (Figure |10[ ) and 
from Schiedea (Figure [ll| . In contrast, the empirical distribution for the highly-inbred BP population of 
K. marmoratus (Figure shows an absence of individuals formed by random outcrossing (T = 0). That 
our method accurately estimates T from simulated data (Figure argues against attributing the inferred 
deficiency of biparental individuals in the BP data set to an artifact of the method. Rather, the deficiency may 
indicate a departure from the underlying model, which assumes reproduction only through self-fertilization 
or random outcrossing. In particular, biparental inbreeding as well as selling may reduce the fraction of 
individuals formed by random outcrossing. Mis-scoring of heterozygotes as homozygotes due to null alleles 
or other factors, a possibility directly addressed by RMES ( David et aZ.|2007 1, may also in principle contribute 
to the paucity of outbred individuals. 


Components of inference 

Locus-specific mutation rates: Our method estimates the scaled mutation rate (|^ at each locus using 
the Dirichlet Process Prior (DPP). This approach improves on existing methods in several ways. First, we 
estimate a single parameter for each locus instead of estimating multiple allele frequencies per locus as do 
Enjalbert and David (20001. Second, we estimate for each locus the scaled mutation rate, a fundamental 
component of the evolutionary process, rather than the heterozygosity 0. a random outcome of that process. 
Third, incorporation of the DPP permits the simultaneous estimation of the number of classes of mutation 
rates, the mutation rate for each class, and the class membership of each locus. It accords the increased 
accuracy derived from pooling loci with similar mutation rates without a priori knowledge of the partitioning 
of loci among rate classes or even the number of classes. 


Joint inference of mutation and inbreeding rates: For the infinite-alleles model of mutation, the Ewens 


Sampling Formula (ESF, Ewens 19721 provides the probability of any allele frequency spectrum (AFS) 


observed at a locus in a sample derived from a panmictic population. Under partial self-fertilization, the 
ESF provides the probability of an AFS observed among genes, each sampled from a distinct individual. 
For such genic (as opposed to genotypic) samples, the coalescence process under inbreeding is identical to 


the standard coalescence process, but with a rescaling of time (Fu 1997 Nordborg and Donnelly 19971. 


Accordingly, genic samples may serve as the basis for the estimation of the single parameter of the ESF, the 
scaled mutation rate 0 * ([^, but not the rate of inbreeding apart from the scaled mutation rate. 

Our method uses the information in a genotypic sample, the genotype frequency spectrum, to infer both 
the uniparental proportion s* and the scaled mutation rate 9*. Our sampler reconstructs the genealogical 
history of a sample of diploid genotypes only to the point of the most recent random-outcross event of each 
individual, with the number of consecutive generations of inbreeding in the immediate ancestry of a given 
individual (T^ for individual k) corresponding to a latent variable in our Bayesian inference framework. 
Invocation of the ESF beyond the point at which all lineages reside in separate individuals obviates the 
necessity of further genealogical reconstruction. As a consequence, our method may be better able to 
accommodate genome-scale magnitudes of observed loci (L). 

Identity disequilibrium ( Cockerham and Weir||[T96^ , the correlation in heterozygosity across loci within 
individuals, reflects that all loci within an individual experience the most recent random-outcross event at the 
same time, irrespective of physical linkage. The heterozygosity profile of individual k provides information 
about Tfe (16), which in turn reflects the uniparental proportion s*. Observation of multiple individuals 


provides a basis for inference of both the uniparental proportion s* and the scaled mutation rate 9*. 


Identifiability: In an analysis based solely on the genotype frequency spectrum observed in a sample, 
the likelihood depends on just two composite parameters: the probability that a random individual is 
uniparental (s*) and the scaled rates of mutation 0 * 
set 


(20) across loci. Any model for which the parameter 


comprises more than one parameter is not fully identifiable from the genetic data alone. In the 
pure hermaphroditism model ([^, for example, basic parameters s (fraction of fertilizations by selling) and 
r (relative viability of uniparental offspring) are nonidentifiable: any assignments that determine the same 
values of composite parameters s* and ©* have the same likelihood. 
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For each basic parameter in ^ beyond one, identifiability requires incorporation of additional information 
beyond the genetic data. A full treatment of such information requires expansion of the likelihood function to 
encompasses an explicit model of the new information. Our androdioecy model (|^, for example, comprises 
3 parameters, including the frequency of males among reproductives {pm) well as s and r. In our analysis 
of microsatellite data from the killifish Kryptolebias marmoratus ( Mackiewicz et ^|2006 Tatarenkov et al. 
20121 , the expanded likelihood function corresponds to the product of the probability of the genetic data and 


the probability of the number of males observed among a total number of individuals (291. In the absence of 


information regarding inbreeding depression (r), we assigned r = 1 to permit estimation of the uniparental 
proportion (s*) under a uniform prior distribution. This assignment does not affect the reliability of our 
estimates (s*, ©*, Pm, Sa, etc.)', rather, the analysis is agnostic concerning the influence of the relative 
viability of inbred offspring (r) and the rate of self-fertilization (s) in determining the probability that a 
random individual is uniparental (s*). 

Non-identifiable parameters can also be estimated through the incorporation of informative priors. Be¬ 
cause identifiability is defined in terms of the likelihood, which is unaffected by priors, such parameters 
remain non-identifiable. Even so, informative priors assist in their estimation through Bayesian approaches. 


which do not require parameters to be identifiable. To explore the data set from Schiedea salicaria (Wallace 


et al. 20111, we use our 4-parameter gynodioecy model (101, the basic parameters of which include the 


proportion of females among reproductives (p/), the relative seed set of females (cr), the relative viability of 
uniparental offspring (t), and the proportion of seeds of hermaphrodites set by self-pollen (a). In a manner 
similar to the androdioecy study, our analysis uses an extended likelihood function, modeling the number 


of females as a binomial random variable (311. In addition, we use earlier experimental evidence to justify 


the assignment of cr = 1 (Weller and Sakai 20051 and to develop an informative prior for r ((30l: Sakai 
et al. 19891. This procedure permits estimation of 3 basic parameters, including the proportion of seeds of 


hermaphrodites set by self-pollen (a). 


Beyond estimation of the selfing rate 

Our MCMC implementation updates the full set of basic parameters, with likelihoods determined from the 
implied values of composite parameters s* and ©*. Incorporation of additional information, either through 
extension of the likelihood or through informative priors, permits inference not only of the basic parameters 
but also of functions of the basic parameters. For example. Table includes estimates of the proportion 
of seeds of hermaphrodites set by self-pollen (a) and the probability that a random gene derives from a 
female parent ((1 — sg)F/2) in gynodioecious S. salicaria. We are not aware of other studies in which these 
quantities have been inferred from the pattern of neutral genetic variation observed in a random sample. 

Among the most biologically-significant functions to which this approach affords access is relative effective 
number S' ([^, a fundamental component of the reproductive value of the sexes ( Fisher|1958 l. We denote the 
probability that a pair of genes, randomly drawn from distinct individuals, derive from the same parent in the 
preceding generation as the rate of parent-sharing (1/A^*). Its inverse (N*) corresponds to the “inbreeding 
effective size” of Crow and Denniston (19881. Relative effective number S is the ratio of N* to the total 


number of reproductive individuals. For example, in the absence of inbreeding (s* = 0), N* in our gynodioecy 
model (10) corresponds to Wright’s (19691 harmonic mean expression for effective population size and S to 
the ratio of N* and Nf + Nh, the total number of reproductive females and hermaphrodites. In general 
(s* > 0), relative effective size S reflects reductions in effective size due to inbreeding in addition to differences 
in numbers of the sexual forms. Figure [represents posterior distributions of S for the 3 data sets explored 
here. These results suggest that relative effective number S in each of the natural populations surveyed 
lies close to its maximum of unity, with the effective number defined through the rate of parent-sharing 
approaching the total number of reproductives. 
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Figure 12 Posterior distributions of relative effective number S Q for data sets derived from Kryptolebias 
marmoratus (BP and TC populations) and Schiedea salicaria. 
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Appendix A The last-sampled gene 


We address the probability that the last-sampled gene in a sample of size i represents a novel allele (22a I. 


Under the infinite alleles model of mutation, a single mutation in a lineage suffices to distinguish a new 
allele. We denote the last-sampled gene in a sample of size i as the focal gene, and consider the level of the 
genealogical tree in which its ancestral lineage either receives a mutation or joins the gene tree of the sample 
at size (i — 1). Level I of the entire (i-gene) gene tree corresponds to the segment in which I lineages persist. 

The probability that the line of descent of the focal gene terminates in a mutation immediately, in level 
i of the genealogy, is 

u _ 0* 

nu+[l)/N* i{9*+i-l)' 

In general, the probability that the lineage of the focal gene terminates on level I > 2 is 


(*-2)u+C-2)/W 


Iu+Q/N* 


lU ■ 


+ {f,)IN* {i-l)u+{^-^)/N* ■■■ (l + l)u+{%^)/N* lu+{]^)/N 


0 * 


i{9* -I- i — 1) 


This expression illustrates the invariance over termination orders noted by Griffiths and Lessard (20051. 


Summing over all levels, including level 2, for which a mutation in either remaining lineage ensures that the 
focal gene represents a novel allele, we obtain the overall probability that the last-sampled gene represents 
a novel allele: 

0 *{i-2) 20* 0* 


i{0* +i — l) i(0* -I- i — 1) 


- i — 1 


Appendix B Estimators of Fjs 


We follow |Weir| ( |1996 1 in developing an estimate of the uniparental proportion s* from Fis alone (281. 

For a single locus, a simple estimator of Fjs corresponds to 

for O the observed fraction of heterozygotes in the sample and E the expected fraction based on Hardy- 
Weinberg proportions given the observed allele frequencies. Explicitly, we have 


Fis = 1 - 






P —if 

J uu Ft 


1 - J2upI 


1 - 


for pu the frequency of allele u in the sample and P^u the frequency of homozygous genotype uu in the 
sample. However, this estimator can be substantially biased for small samples, leading to underestimation 
of Fis ( Weir|[T996 |. 

To address this bias and accommodate multiple loci, we instead adopt 


Fis = 


Eti [EuU {Piuu-Pl) + (l - EuU Piuu) /2n 


Ef=i [(l - EuUpI) - (l - EuLi Piuu) /2r 


(B.l) 


for n the number of diploid genotypes observed, L the number of loci, and Ki the number of alleles at 
locus 1. While this estimator is also biased in general, it corresponds to the ratio of unbiased estimators of 
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FiS ■ El (1 -E,..pI,) a nd (l " EuPfu)^ in which piu is the frequency of allele u at locus I in the entire 
population ( |Weir|[l9^ I. Our analysis of simulated data (Appendix indicates that this estimator is more 
accurate than an estimator that simply averages single-locus estimates: 


- 1 ^ EuLi {Piuu Pi^ + EuLi Piuu^ 


1^1 


{^-Et.Pl) - {l-EuUPiuu) /2n 


(B.2) 


Our A/s-based estimates (281 incorporate (B.ll and not (B.2). 


Appendix C Implementation of the MCMC 


State space: The state space for the Markov chain of our MCMC sampler includes times across sampled 


individuals since the last outcross event T (14), coalescence events across individuals and loci since that event 


I (151, and model-specific parameters (24). The state space also comprises the scaled mutation rates ©* 


(201, which are determined by C, a list specifying the mutation rate category Ci for locus I = 1... L, and Z, 


a list specifying the scaled mutation rate Zi for category f = 1... L -|- 4. In particular, the scaled mutation 
rate at locus I corresponds to 

o: = zc,. (c.i) 

At any given point in the MCMC, the state of the Markov chain corresponds to (I, T, \I/, C, Z). 


Iterations: Each iteration of our MCMC sampler performs multiple updates, with each variable updated 
at least once per iteration. We recorded the state sampled by the MCMC at each iteration. For analyses of 
simulated data sets, we ran Markov chains for 2000 iterations, discarding the first 200 iterations as burn-in. 
For analyses of the actual data sets, we ran Markov chains for 100,000 iterations, discarding the first 10,000 
iterations as burn-in. Convergence appeared to occur as rapidly for actual data as for simulated data, but 
we found empirically that the larger number of samples were needed to achieve smooth density plots for the 
actual data sets. 


Transition kernels: Updating of the continuous variables of mutation rates {Zi} (C.I I and model-specific 
parameters (241 uses both Metropolis-Hastings (MH) transition kernels and auto-tuned slice-sampling 
transition kernels. Updating of the discrete variables {Ci} uses a Gibbs transition kernel. 


Efficient inference on selfing times through collapsed Metropolis-Hastings: Simple Metropolis-Hastings 
(MH) proposals that separately update the time since the most recent outcross event (T^) and coalescence 
history since that event (/.fc) lead to extremely poor mixing efficiency. Strong correlations between and 
I.k cause changes to to be rejected with high probability unless I.k is updated as well. For example, 
consider proposing a change of Tk from 1 to 0. When = 1, on average Iik will be 1 at half of the loci and 0 
at the remaining loci. If any of the Iik = 1, a move to = 0 will always be rejected because the probability 
of a coalescence event more recently than the most recent outcross event is 0 if the sampled individual is 
itself a product of outcrossing. To permit acceptance of changes to Tk, we introduce a proposal for Tk that 
also changes I.k- 

The scheme starts from the value Tk = tk and proposes a new value tj.. In standard MH within Gibbs, 
we would compute the probability of Tk = tk and of Tk = t'f, given that all other parameters are unchanged. 
We modify this MH scheme to compute probabilities without conditioning on the coalescence indicators for 
individual k. However, the coalescence indicators for other individuals are still held constant. To compute 
this probability, let J indicate all the coalescence indicators I.y where y ^ k. Then 

Pr(X,T,J,s,6l) = Pr(X,J|T,s,6»)Pr(T|s)Pr(s)Pr(6l). 

We introduce 1.^ by summing over all possible values i.^. 

Pr(X,J|T,s,0) = ^Pr(X,I.fc =i.fc,J|T,s,0). 

^-fc 
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Since the iik for different loci are independent given T^, we have 

L 

Pr(X,J|T,s,0) = ^[]Pr(X,,/,fe = f,fc,J,|T,s,0) 

^ • fc ^— 1 

L 

1 — 1 ilk 

Therefore, for specific values of T and J, we can compute the sum over all possible values of I.^ for Z = 1... L 
in computation time proportional to L instead of 2^. This is possible because the L coalescence indicators 
for individual k each affect different loci, and are conditionally independent given Tk and J. 

After accepting or rejecting the new value of with I.k integrated out, we must choose new values for 
I.fc given the chosen value of T^. Because of their conditional independence, we may separately sample each 
coalescence indicator Afc for Z = 1... L from its full conditional given the chosen value of Tfc. This completes 
the collapsed MH proposal. 


Appendix D Analysis of simulated data 

Simulations: Our simulator (https://github.com/skumagai/selfingsim) was developed using simuPOP, 
publicly available at http://simupop.sourceforge.net/. It explicitly represents N = 10,000 individuals, 
each bearing two genes at each of L unlinked loci. Mutations arise at locus Z at scaled rate 9i 0 , in 
accordance with the the infinite-alleles model. 

We assigned to uniparental proportion s* values ranging from 0.01 to 0.99, with half of the L = 32 loci 
assigned scaled mutation rate 6 = 0.5 and the remaining loci 9 = 1.5. 

We conducted 10^ independent simulations for each assignment of s*. Each simulation was initialized 
with each of the 27V x 32 genes representing a unique allele. Most of this maximal heterozygosity was lost very 
rapidly, with allele number and allele frequency spectrum typically stabilizing well within 107V generations. 
After 207V generations, we recorded the realized population, from which 100 independent samples of L = 32 
loci of size n = 70 were extracted. From this collection, we randomly chose L = 6 loci and subsampled 100 
independent samples of size n = 6. 


Analysis: To 10^ independent samples from each of 10^ independent simulations for each assignment of 
the uniparental proportion s*, we applied our Bayesian method, the Fjs method, and RMES. Our Bayesian 
method is open-source and can be obtained at 


https://github.com/bredelings/BayesianEstimatorSelfing/ 


We used the implementation of RMES (]David et aZ.||2007 1 provided at 


http://www.cefe.cnrs.fr/images/stories/DPTEEvolution/Genetique/fichiers7o20Equipe/RMES7, 
2020097.2827.29.zip 
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Supplementary Methods 


1 Indicators of accuracy 

To compare the accuracy of our Bayesian method to RMES and the Fjs method, which prodnce point esti¬ 
mates, we summarize the posterior distribution of the uniparental proportion s* by the median. Here, we 
compare the median to the mode and mean of the posterior distribution. 

Figure [^suggests that the bias and root-mean-squared (rms) error of these three indices exhibit different 
properties. For example, the posterior mode shows smaller bias thronghont the parameter range, bnt the 
median and mean show smaller rms error for s* near the bonndaries (near 0 or 1). 


2 Average error 


As for the case of large simulated data sets (Figure]^, Figure S2 indicates that upon application to smaller 
samples (n = 10 individuals, L — Q loci), both RMES and onr method show positive bias npon application to 
data sets for which the true uniparental proportion s* is close to zero and negative bias for s* close to unity. 
It further indicates that while both methods exhibit more error for small samples than large samples, onr 
Bayesian method exhibits less error than RMES thronghont the range of the nniparental proportion (s*). 


3 Frequent ist coverage 


As for the 95% BCIs (Fignrej^, Fignre S3 indicates that BCIs of different nominal 
0.95, and 0.99) display the same pattern, with coverage exceeding the desired value 
s* values and dipping below the desired value for very high values of s*. Coverage is 
value for the 0.99 and 0.95 levels than for the 0.5 level. 


valnes (0.5, 0.75, 0.9, 
for intermediate trne 
closer to the nominal 


4 Data analysis 


4.1 Androdioecious vertebrate 

Low outcrossing rate: As noted in the main text, we find evidence of a multimodal distribution of mutation 
rates in the BP population of K. marmoratus. 

Figure [S5| shows the posterior distributions of number of generations since the most recent outcross event 

(Hi- 


Higher outcrossing rate: Figure S6 presents posterior distributions of locus-specific mutation rates (com¬ 
pare Fignre S4|. For each individual in the TC sample. Figure S7 shows the posterior distribution of the 
number of consecutive generations of selling in its immediate ancestry. 


4.2 Gynodioecious plant 

Figure |58] presents posterior distributions for locus-specific mutation rates inferred from the S. salicaria data 
set. The loci appear to have similar posterior medians. 

Figure |S9| presents the inferred number of generations since the most recent outcross event (141 for 
each individual k. 

Figure SIO presents posterior distributions for the uniparental proportion (sg), the proportion of females 
among reprodnctives (p/), the proportion of seeds set by hermaphrodites by self-pollen (a), and the viability 
of uniparental offspring relative to biparental offspring (r). 
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Figure SI Errors for the posterior mean, posterior median, and posterior mode. Blue curves (rms) indicate 
the root-mean-squared error, and red curves (bias) the average deviation. Averages are taken across simulated 
data sets at each true value of the selling rate s*. 
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Figure S2 Errors for the full likelihood (posterior median), RMES, and Fjs methods for a small sample 
(n = 10 individuals, L = 6 loci). In the legend, rms indicates the root-mean-squared error and bias the 
average deviation. Averages are taken across simulated data sets at each true value of s*. 
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Figure S3 Frequentist coverage for Bayesian credible intervals at different levels of credibility under the large 
sampling regime (n = 70, L = 32). 
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Figure S4 Posterior distributions for mutation rates at each locus in K. marmoratus (BP population). For 
each distribution, the locus name is indicated in the grey shaded box. 
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Figure S5 Number of generations since the most recent outcross event in the ancestry of each individual in the 
sample from the BP population of K. marmoratus. The area of each dot indicates the posterior probability 
that an individual (X-axis) has the indicated number (Y-axis) of consecutive generations of selling in its 
immediate ancestry. The blue line indicates the posterior mean number of selling generations and the red 
line indicates the number of heterozygous loci across individuals. The Y-axis is truncated to [0,30]. 
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Figure S6 Mutation rates at each locus for K. marmoratus (TC population). For each distribution, the 
locus name is indicated in the grey shaded box. 
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Figure S7 Number of generations since the most recent outcross event in the ancestry of each individual in 


the sample from the TC population of K. marmoratus. Symbols as in Figure S5 
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Figure S8 Posterior distributions for mutation rates at locns in S. salicaria. For each distribution, the locns 
name is indicated in the grey shaded box. 
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Figure S9 Estimated number of selfing generations for each individual for S. salicaria. The area of each dot 
indicates the posterior probability that a numbered individual (x-axis) has been selfed for a given number 
of generations (y-axis). For each individual the blue line indicates the posterior mean number of selfing 
generations and the red line indicates the nnmber of heterozygons loci. 
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Figure SIO Posterior distributions on (a) sg, (b) pf, (c) a, and (d) r for the Schiedea salicaria data set. 
Also shown are 95% BCI (red), 50% BCI (orange), and median (black dot). 
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